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The  purpose  of  this  article  is  numerical  verification  of  the  theory  of  weak  turbulence. 
We  performed  numerical  simulation  of  an  ensemble  of  nonlinearly  interacting  free  grav¬ 
ity  waves  (swell)  by  two  different  methods:  solution  of  primordial  dynamical  equations 
describing  potential  flow  of  the  ideal  fluid  with  a  free  surface  and,  solution  of  the  kinetic 
Hasselmann  equation,  describing  the  wave  ensemble  in  the  framework  of  the  theory  of 
weak  turbulence.  In  both  cases  we  observed  effects  predicted  by  this  theory:  frequency 
downshift,  angular  spreading  and  formation  of  Zakharov-Filonenko  spectrum  ~ 

To  achieve  quantitative  coincidence  of  the  results  obtained  by  different  methods,  one  has 
to  supply  the  Hasselmann  kinetic  equation  by  an  empirical  dissipation  term  S^iss  modeling 
the  coherent  effects  of  white-capping.  Using  of  the  standard  dissipation  terms  from  op¬ 
erational  wave  predicting  model  ( WAM)  leads  to  significant  improvement  on  short  times, 
but  not  resolve  the  discrepancy  completely,  leaving  the  question  about  optimal  choice  of 
Sdiss  open.  In  a  long  run  WAM  dissipative  terms  overestimate  dissipation  essentially. 

1.  Introduction. 

The  theory  of  weak  turbulence  is  designed  for  statistical  description  of  weakly-nonlinear 
wave  ensembles  in  dispersive  media.  The  main  tool  of  weak  turbulence  theory  is  kinetic 
equation  for  squared  wave  amplitudes,  or  a  system  of  such  equations.  Since  the  discovery 
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of  the  kinetic  equation  for  bosons  by  Nordheim  [[T]  (see  also  paper  by  Peierls  [[2])  in  the 
context  of  solid  state  physics,  this  quantum-mechanical  tool  was  applied  to  wide  variety  of 
classical  problems,  including  wave  turbulence  in  hydrodynamics,  plasmas,  liquid  helium, 
nonlinear  optics,  etc.  (see  monograph  by  Zakharov,  Falkovich  and  LVov  [|3]).  Such  kinetic 
equations  have  rich  families  of  exact  solutions  describing  weak-turbulent  Kolmogorov 
spectra.  Also,  kinetic  equations  for  waves  have  self-similar  solutions  describing  temporal 
or  spatial  evolution  of  weak  -  turbulent  spectra. 

However,  in  our  opinion,  the  most  remarkable  example  of  weak  turbulence  is  wind- 
driven  sea.  The  kinetic  equation  describing  statistically  the  gravity  waves  on  the  surface 
of  ideal  liquid  was  derived  by  Hasselmann  [  H] .  Since  this  time  the  Hasselmann  equation 
is  widely  used  in  physical  oceanography  as  foundation  for  development  of  wave-prediction 
models  such  as  WAM,  SWAN  and  WAVEWATCH:  it  is  quite  special  case  between  other 
applications  of  the  theory  of  weak  turbulence  due  to  the  strength  of  industrial  impact. 

In  spite  of  tremendous  popularity  of  the  Hasselmann  equation,  its  validity  and  appli¬ 
cability  for  description  of  real  wind-driven  sea  has  never  been  completely  proven.  It  was 
criticized  by  many  respected  authors,  not  only  in  the  context  of  oceanography.  There  are 
at  least  two  reasons  why  the  weak-turbulent  theory  could  fail,  or  at  least  be  incomplete. 

The  hrst  reason  is  presence  of  the  coherent  structures.  The  weak-turbulent  theory 
describes  only  weakly-nonlinear  resonant  processes.  Such  processes  are  spatially  extended, 
they  provide  weak  phase  and  amplitude  correlation  on  the  distances  significantly  exceeding 
the  wave  length.  However,  nonlinearity  also  causes  another  phenomena,  much  stronger 
localized  in  space.  Such  phenomena  -  solitons,  quasi-solitons  and  wave  collapses  are 
strongly  nonlinear  and  cannot  be  described  by  the  kinetic  equations.  Meanwhile,  they 
could  compete  with  weakly-nonlinear  resonant  processes  and  make  comparable  or  even 
dominating  contribution  in  the  energy,  momentum  and  wave-action  balance.  For  gravity 
waves  on  the  fluid  surface  the  most  important  coherent  structures  are  white-cappings  (or 
wave-breakings),  responsible  for  essential  dissipation  of  wave  energy. 

The  second  reason  limiting  the  applicability  of  the  weak-turbulent  theory  is  finite  size 
of  any  real  physical  system.  The  kinetic  equations  are  derived  only  for  infinite  media, 
where  the  wave  vector  runs  continuous  d-dimensional  Fourier  space.  Situation  is  different 
for  the  wave  systems  with  boundaries,  e.g.  boxes  with  periodical  or  reflective  boundary 
conditions.  The  Fourier  space  of  such  systems  is  infinite  lattice  of  discrete  eigen-modes. 
If  the  spacing  of  the  lattice  is  not  small  enough,  or  the  level  of  Fourier  modes  is  not  big 
enough,  the  whole  physics  of  nonlinear  interaction  becomes  completely  different  from  the 
continuous  case.  We  shall  call  effects  caused  by  a  finite  size  of  a  system  “mesoscopic 
effects”.  These  effects  could  be  important  in  nature  and  they  certainly  should  be  taken 
in  to  account  by  anyone  who  performs  numerical  simulation  of  wave  turbulence. 

For  these  two  reasons  verihcation  of  the  weak  turbulent  theory  is  an  urgent  problem, 
important  for  the  whole  physics  of  nonlinear  waves.  The  verihcation  can  be  done  by  direct 
numerical  simulation  of  the  primordial  dynamical  equations  describing  wave  turbulence 
in  nonlinear  medium. 

So  far,  the  numerical  experimentalists  tried  to  check  some  predictions  of  the  weak- 
turbulent  theory,  such  as  weak-turbulent  Kolmogorov  spectra.  For  the  gravity  wave  tur¬ 
bulence  the  most  important  is  Zakharov-Filonenko  spectrum  ~  w--*  [  a.  At  the 
moment,  this  spectrum  was  observed  in  numerous  numerical  experiments  [IB]-[I^. 
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The  attempts  of  verification  of  weak  turbulent  theory  through  numerical  simulation  of 
primordial  dynamical  equations  has  been  started  with  numerical  simulation  of  2D  optical 
turbulence  [EH,  which  demonstrated,  in  particular,  co-existence  of  weak-turbulent  and 
coherent  events. 

Numerical  simulation  of  2D-turbulence  of  capillary  waves  was  done  in  10.  lEl.  and 
[H].  The  main  results  of  the  simulation  consisted  in  observation  of  classical  regime  of 
weak  turbulence  with  spectrum  ~  and  discovery  of  non-classical  regime  of 

“frozen  turbulence”,  characterized  by  absence  of  energy  transfer  from  low  to  high  wave- 
numbers.  The  classical  regime  of  turbulence  was  observed  on  the  grid  of  256  x  256  points 
at  relatively  high  levels  of  excitation,  while  the  “frozen”  regime  was  realized  at  lower 
levels  of  excitation,  or  more  coarse  grids.  The  effect  of  “frozen”  turbulence  is  explained 
by  mesoscopic  effects:  sparsity  of  both  exact  and  approximate  resonances.  The  classical 
regime  of  turbulence  becomes  possible  due  to  broadening  of  resonances  by  nonlinearity. 
The  classical  and  the  frozen  regime  could  coexist.  We  call  this  situation  “mesoscopic 
turbulence” . 

In  fact,  the  “frozen”  turbulence  is  close  to  KAM  regime,  when  the  dynamics  of  turbu¬ 
lence  is  close  to  the  behavior  of  integrable  system  [  [8] . 

Simulation  of  the  surface  gravity  waves  turbulence  for  the  hrst  time  was  done  simultane¬ 
ously  by  Tanaka  m  and  Onorato  at  ah  [HQ].  Due  to  limitations  of  calculation  performance 
of  that  time  computers  simulations  were  limited  to  dynamical  equations  and  quite  short 
simulation  times.  After  that  sime  many  articles  developed  these  pioneering  results  [m 
HH.  The  nonlinear  effects  for  gravity  waves  are  weaker  than  for  cappilary  waves,  as  a 
result  the  influence  of  the  mesoscopic  effects  is  stronger.  It  makes  the  simulations  much 
more  difficult.  In  the  experiments  mentioned  above,  the  grid  was  hne  enough  to  resolve 
the  spectral  tails  and  observe  the  weak  turbulent  Kolmogorov  asymptotics  Ik  ~  k~^,  but 
it  was  too  coarse  to  avoid  mesoscopic  effects  in  the  area  of  spectral  peak.  The  hnest  grid 
2048  X  4096  was  used  by  Tanaka,  however  in  his  experiments  the  spectral  peak  was  posed 
at  k  =  68,  while  the  level  of  nonlinearity  was  quite  small  (typical  steepness  ju  ~  0.07). 
More  over,  his  experiment  was  very  short  in  terms  of  characteristic  time  of  the  waves  (only 
several  tens  periods  of  the  leading  wave).  If  he  would  continue  his  calculations  he  would 
observe  formation  of  the  weak  turbulent  Kolmogorov  tail  as  in  the  work  by  Onorato  at 
al.  [  [To] .  However  his  spectral  peak  was  seriously  damaged  by  mesoscopic  effects. 

In  this  article  we  present  the  results  of  simulation  of  the  surface  gravity  waves  turbu¬ 
lence,  namely  modelling  of  swell  propagation.  All  previous  experimentalists  since  Tanaka  [ 
E]  and  Onorato  at  al.  [□TO]  performed  only  one  type  of  experiments  —  solution  of  the  pri¬ 
mordial  dynamical  equations.  We  made  two  experiments  simultaneously.  We  solved  not 
only  dynamical  equations  but  also  Hasselmann  kinetic  equation  and  compared  the  results. 
We  think  that  our  results  can  be  considered  as  the  hrst  attempt  of  direct  verihcation  of 
Hasselmann  kinetic  equation. 

In  the  dynamic  experiment  we  used  less  hne  grid  than  Tanaka  (512  x  4096)  but  we  posed 
the  spectral  peak  at  k  =  300,  and  our  initial  steepness  was  much  higher  than  in  the  Tanaka 
case  {fi  ~  0.15).  Due  to  much  higher  k  and  stronger  nonlinearity  we  managed  basically  to 
suppress  the  mesoscopic  ehects.  The  bulk  of  energy  containing  modes  satishes  Rayleigh 
distribution,  according  to  the  weak-turbulent  scenario.  The  distribution  has  a  heavy  tail 
of  abnormally  intensive  harmonics  (“oligarchs”  according  to  terminology,  introduced  in 
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the  article  mg).  but  they  contain  no  more  than  5%  of  total  energy. 

One  important  point  should  be  mentioned.  In  our  experiments  we  observed  not  only 
weak  turbulence,  but  also  additional  nonlinear  dissipation  of  the  wave  energy,  which  could 
be  identihed  as  the  dissipation  due  to  white-capping.  We  observed  fast  broadening  of  the 
spectra  due  to  multiple  harmonics  generation.  Formation  of  the  broad  spectrum  with 
strong  small  scale  tails  can  be  explained  and  formation  of  quite  sharp  structures  on  waves 
crests.  Further  development  of  these  tails  is  suppressed  by  an  artihcial  dissipation  used  in 
our  dynamical  experiment.  One  can  say  that  in  such  a  way  we  roughly  simulated  white¬ 
capping  phenomenon  (continuous  in  time  dissipation  of  energy  due  to  multiple  acts  of 
wave  braking  on  the  very  edge  of  wave  crest,  arresting  formation  of  derivative  singularity 
on  the  sharp  wave  crest).  Of  course  we  cannot  simulate  directly  wave  braking  phenomenon 
but  one  can  say  that  our  method  provide  may  be  rough  but  reasonable  simulation  of  this 
effect. 

To  reach  an  agreement  with  dynamic  experiments,  we  had  to  add  to  the  kinetic  equation 
a  phenomenological  dissipation  term  Sdiss-  In  this  article  we  examined  dissipation  terms 
used  in  the  operational  wave-prediction  models  WAM  Cycle  3  and  WAM  cycle  4  (hereafter 
referenced  as  WAM3  and  WAM4  correspondingly).  Both  of  these  terms  overestimate 
nonlinear  dissipation  signihcantly.  Term  given  in  WAM3  gives  acceptable  results  on  short 
periods  of  time  (time  less  than  IOOOTq,  where  Tq  is  the  time  period  of  the  leading  wave 
of  the  initial  condition).  But  at  the  end  of  simulation  {t  =  3378To)  an  error  in  wave 
action  become  close  to  30%.  For  WAM4  the  situation  is  even  worse.  If  the  characteristic 
wave  length  of  initial  conditions  is  equal  to  22m,  then  3378To  is  a  little  bit  more  than 
only  3  hours  of  wave  development.  Even  primitive  viscous  dissipative  term  without  any 
simulation  of  strongly  nonlinear  events  gives  us  better  results.  It  means  that  the  question 
about  a  reasonable  formulae  for  Sdiss  is  open  for  now. 

2.  Deterministic  and  statistic  models. 

In  the  ’’dynamical”  part  of  our  experiment  surface  of  the  fluid  was  described  by  two 
functions  of  horizontal  variables  x,y  and  time  t:  surface  elevation  r]{x,y,t)  and  velocity 
potential  on  the  surface  'il^i^x,  y,  t).  They  satisfy  the  canonical  equations  [  l24] 

dr]  6H  d'lp  6H 

dt  d'lp'  dt  dr]' 

Hamiltonian  H  is  presented  by  the  hrst  three  terms  in  expansion  on  powers  of  nonlinearity 
V?7 

H  =  H0  +  H1  +  H2  + 


f 


n  I'lV'i/'P  —  dxdv. 


(2) 


Here  k  is  the  linear  integral  operator  k  =  V— V^,  dehned  in  Fourier  space  as 


(3) 
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Using  Hamiltonian  ([2])  and  equations  ([T])  one  can  get  the  dynamical  equations  [E]: 


V  = 


k'lp  —  (V  {riV 'ip))  —  k[r]k'ip]  + 
+k{rik[rikip])  +  |V^[r7^fc-^]  + 


^k[T]^V^7p]  +  F-^[^kVk], 

- 


-9V  j  ^ 
-[kip]k[r]kip] 


(4) 


Here  F~^  corresponds  to  inverse  Fourier  transform.  We  introduced  artificial  dissipative 
terms  F~^['ykVk\  and  F~^['jk'ipk],  corresponding  to  pseudo- viscous  high  frequency  damping 
following  recent  work  [  [19] . 

The  model  (II])- (01)  was  used  in  the  numerical  experiments  [E]-[E],  [112],  [03],  [05],  [ 

07],  [0H]. 

Introduction  of  the  complex  normal  variables 


(5) 


where  Uk  =  a/^,  transforms  equations  ([T])  into 

ggg  _  .6H 

dt 

k 


(6) 


To  proceed  with  statistical  description  of  the  wave  ensemble,  hrst,  one  should  perform 
the  canonical  transformation  >  6^,  which  excludes  the  cubical  terms  in  the  Hamilto¬ 
nian.  The  details  of  this  transformation  can  be  found  in  the  paper  by  Zakharov  (1999)  [ 
1^.  After  the  transformation  the  Hamiltonian  takes  the  form 


H  =  J  +\J  ^kk^hWAh,  X  (7) 

x%+fc^_fc^_fc3dfcidfc2d4- 

where  T  is  a  homogeneous  function  of  the  third  order: 

T{ek,  eki,  eh,  eh)  =  e^T{k,  h,  h,  h).  (8) 


Connection  between  and  together  with  explicit  expression  for  can  be  found, 

for  example,  in  [125]. 

Let  us  introduce  the  pair  correlation  function 

<  >=  gN^6{k  -  k'),  (9) 

where  is  the  spectral  density  of  the  wave  function.  This  dehnition  of  the  wave  action 
is  common  in  oceanography. 

We  also  introduce  the  correlation  function  for  transformed  normal  variables 

<  gnfS{k-k') 


(10) 
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Functions  and  can  be  expressed  through  each  other  in  terms  of  cumbersome  power 
series  [  |25]  of  expansion  on  p.  Here  fi  is  the  characteristic  steepness  dehned  as  follows 

where  E  is  the  wave  energy  and  N  is  the  wave  action.  Following  this  dehnition  for  the 
Stokes  wave  of  small  amplitude 


Tj  =  a  cos(/cx), 
fi  ~  ak. 

On  deep  water  their  relative  difference  between  and  Nj^  is  of  the  order  of  and  can 
be  neglected  (in  most  cases  experimental  results  shows  /i  ~  0.1). 

Spectrum  satishes  Hasselmann  (kinetic)  equation  [  H] 

dTi~* 

*5*77,/ [ti]  “1“  Sdiss  “1“ 

Sni[n]  =  2TTg‘^  j  + 

\  (12) 
-  ^kVk2  -  ^kVks)  X 
xS  {Uk  +  UJki  —  tOk2  ~  ^ki)  X 
x5  {k  +  ki  —  k2  —  ksj  dkidk2dks. 

Here  Sdiss  is  an  empiric  dissipative  term,  corresponding  to  white-capping. 

Stationary  conservative  kinetic  equation 


^n;  =  0 


(13) 


has  the  rich  family  of  Kolmogorov-type  [  [26]  exact  solutions.  Among  them  is  Zakharov- 
Filonenko  spectrum  [[5]  for  the  direct  cascade  of  energy 

Uk  ~  (14) 

and  Zakharov- Zaslavsky  [ET],  [EH]  spectra  for  the  inverse  cascade  of  wave  action 


nk  ~ 


1 

fc23/6  ’ 


(15) 


3.  Deterministic  Numerical  Experiment. 

3.1.  Problem  Setup 

The  dynamical  equations  (jl])  have  been  solved  in  the  real-space  domain  27r  x  27r  on 
the  grid  512  x  4096  with  the  gravity  acceleration  set  to  =  1.  The  solution  has  been 
performed  by  the  spectral  code,  developed  in  [E2]  and  previously  used  in  [E3],[[l2],  [113]. [ 
[15].  We  have  to  stress  out  that  in  the  current  computations  the  resolution  in  K-direction 
(long  axis)  is  better  than  the  resolution  in  X-direction  by  the  factor  of  8. 

This  approach  is  reasonable  if  the  swell  is  essentially  anisotropic,  almost  one-dimensional. 
This  assumption  will  be  validated  by  the  proper  choice  of  the  initial  data  for  computation. 
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As  the  initial  condition,  we  used  the  Gaussian-shaped  distribution  in  Fourier  space  (see 
Fig.  [I]): 


I  Or  I  =  Ajexp 


1  k  —  kc 

2  DJ 


m  =  10 


-12 


V 

k  —  kr 


k  —  kr 


<2A, 


>2A, 


Ai  =  0.92  X  10"®,  A  =  60, 
ko  =  (0;300),a;o  = 


(16) 


The  initial  phases  of  all  harmonics  were  random.  The  average  steepness  of  this  initial 
condition  was  /r  ~  0.15  (dehned  in  accordance  with  EqJTTll. 

To  realize  similar  experiment  in  the  laboratory  wave  tank,  one  has  to  to  generate  the 
waves  with  wave-length  300  times  less  than  the  length  of  the  tank.  The  width  of  the  tank 
would  not  be  less  than  1/8  of  its  length.  The  minimal  wave  length  of  the  gravitational  wave 
in  absence  of  capillary  effects  can  be  estimated  as  A^m  —  3cm.  The  leading  wavelength 
should  be  higher  by  the  order  of  magnitude  A  ~  30cm. 

In  such  big  tank  of  200  x  25  meters  experiment ators  can  observe  the  evolution  of  the 
swell  until  approximately  700To  -  still  less  than  in  our  experiments.  In  the  tanks  of  smaller 
size,  the  effects  of  discreetness  the  Fourier  space  will  be  dominating,  and  experimentalists 
will  observe  either  “frozen”,  or  “mesoscopic”  wave  turbulence,  qualitatively  different  from 
the  wave  turbulence  in  the  ocean. 

To  stabilize  high-frequency  numerical  instability,  the  damping  function  has  been  chosen 
as 


^  j  0,k  <  kd, 

\  -j{k-kdf,k>kd,  (17) 

kd  =  1024,7  =  5.65  x  10"^ 

The  simulation  was  performed  until  t  =  1225,  which  is  equivalent  to  3378  Tq,  where 
To  =  27r/-y/A  is  the  period  of  the  wave,  corresponding  to  the  maximum  of  the  initial 
spectral  distribution. 


3.2.  Zakharov-Filonenko  spectra 

Like  in  the  previous  papers  [  [iQ] ,  [  [l2] ,  [  |T3]  and  [HH],  we  observed  fast  formation  of  the 
spectral  tail,  described  by  Zakharov-Filonenko  law  for  the  direct  cascade  ^  k  ^  m 
(see  Figj2]).  In  the  agreement  with  [[I5],  the  spectral  maximum  slowly  down-shifts  to  the 
large  scales  region,  which  corresponds  to  the  inverse  cascade  [I?F].[I28]. 

Also,  the  direct  measurement  of  energy  spectrum  has  been  performed  during  the  final 
stage  of  the  simulation,  when  the  spectral  down  shift  was  slow  enough.  This  experiment 
can  be  interpreted  as  the  ocean  buoy  record  -  the  time  series  of  the  surface  elevations  has 
been  recorded  at  one  point  of  the  surface  during  Tbuoy  —  300To.  The  Fourier  transform  of 
the  autocorrelation  function 


Tbuoy  1'^ 


-Tbuoy!’^ 


<  i](t  +  T)y(T)  > 


(18) 
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Figure  1.  Initial  distribution  of 


on  fc-plane. 


allows  to  detect  the  direct  cascade  spectrum  tail  proportional  to  (see  Figj3]),  well 
known  from  experimental  observations  [  [29] ,  [  |30] ,  [  |3T] . 

In  this  paper  we  had  no  intention  to  improve  results  of  the  previous  papers.  Inertial 
interval  of  the  angle-averaged  spectra  (cu-spectrum  also  is  angle-averaged  because  it  de¬ 
pends  only  on  frequency  oo  ~  y/k)  is  limited  due  to  anisotropy  of  the  integration  domain. 
In  this  case  we  have  less  than  one  decade  interval  and  we  cannot  hnd  an  exponent  of  the 
Kolmogorov  tail  surely,  however  this  work  was  already  done  in  several  previous  articles 
(see  for  instance  [US]).  Here  we  limit  our  observations  by  qualitative  correspondance. 

3.3.  Is  the  weak-turbulent  scenario  realized? 

Presence  of  Kolmogorov  asymptotic  in  spectral  tails,  however,  is  not  enough  to  validate 
applicability  of  the  weak-turbulent  scenario  for  description  of  wave  ensemble.  We  have  also 
be  sure  that  statistical  properties  of  the  ensemble  correspond  to  weak-turbulent  theory 
assumptions. 

One  should  stress  out  that  at  the  beginning  of  our  experiments  is  a  smooth  function 
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Figure  2.  Angle-averaged  spectrum  =<  >  in  a  double  logarithmic  scale.  The  tail 

of  distribution  hts  to  Zakharov-Filonenko  spectrum. 


of  k.  Only  the  phases  of  the  individual  waves  are  random.  As  shows  numerical  simulation, 
the  initial  condition  flTbl)  (see  Fig{T])  does  not  preserve  its  smoothness  -  it  becomes  rough 
virtually  immediately  (see  FigH]).  The  picture  of  this  roughness  is  remarkably  preserved 
in  many  details,  even  as  the  spectrum  down-shifts  as  a  whole.  This  roughness  does 
not  contradict  the  weak-turbulent  theory.  According  to  this  theory,  the  wave  ensemble 
is  almost  Gaussian,  and  both  real  and  imaginary  parts  of  each  separate  harmonics  are 
not-correlated.  However,  also  according  to  the  weak-turbulent  theory,  the  spectra  must 
become  smooth  after  averaging  over  long  enough  time  of  more  than  l//i^  periods.  Earlier 
we  observed  such  restoration  of  the  smoothness  in  the  numerical  experiments  of  the  MMT 
model  (see  [I1Z],[I1H],  [09]  and  [|50]).  However,  in  the  experiments  discussed  in  the  article, 
the  roughness  still  persists  and  the  averaging  does  not  suppresses  it  completely.  It  can  be 
explained  by  sparsity  of  the  resonances. 

Resonant  conditions  are  dehned  by  the  system  of  equations: 

T  ojk^  =  ujk2  T 

k  +  k^  =  k2  +  h,  ^  ’ 

These  resonant  conditions  dehne  hve-dimensional  hyper-surface  in  six-dimensional  space 
k,  ki,  /c2-  In  a  hnite  system,  (11911  turns  into  Diophantine  equation.  Some  solutions  of  this 
equation  are  known  [132],  [[n].  In  reality,  however,  the  energy  transport  is  realized  not  by 
exact,  but  approximate  resonances,  posed  in  a  layer  near  the  resonant  surface  and  dehned 
by 

\^k  +  ^ki  —  ^k2  ~  ^k+ki-k2  I  <  T, 


(20) 
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CO 


Figure  3.  Energy  spectrum  in  a  double  logarithmic  scale.  The  tail  of  distribution  fits  to 
asymptotic 


where  F  is  a  characteristic  inverse  time  of  nonlinear  interaction. 

In  the  finite  systems  k,  ki,  ^2  take  values  on  the  nodes  of  the  discrete  grid.  The  weak 
turbulent  approach  is  valid,  if  the  density  of  discrete  approximate  resonances  inside  the 
layer  fl20|)  is  high  enough.  In  our  case  the  lattice  constant  is  Ak  =  1,  and  typical  relative 
deviation  from  the  resonance  surface 


Ao; 

UJ 


00 


oo 


1 


~  2  X  10 


(21) 


Inverse  time  of  the  interaction  F  can  be  estimated  from  our  numerical  experiments:  wave 
amplitudes  change  essentially  during  30  periods,  and  one  can  assume:  F/cn  ~  10“^  S> 

It  means  that  the  condition  for  the  applicability  of  weak  turbulent  theory  is  typically 
satisfied,  but  the  reserve  for  their  validity  is  rather  modest.  As  a  result,  some  particular 
harmonics,  posed  in  certain  “privileged”  point  of  fc-plane  could  form  a  network  of  almost 
resonant  quadruplets  and  realize  signihcant  part  of  energy  transport.  Amplitudes  of  these 
harmonics  exceed  the  average  level  essentially.  This  effect  was  described  in  the  article  [ 
115] ,  where  such  “selected  few”  harmonics  were  called  “oligarchs” .  If  oligarchs  realize  most 
part  of  the  energy  flux,  the  turbulence  is  mesoscopic,  not  weak. 


3.4.  Statistics  of  the  harmonics 

According  to  the  weak-turbulent  scenario,  statistics  of  the  in  any  given  k  should 
be  close  to  Gaussian.  It  presumes  that  the  PDF  for  the  squared  amplitudes  is 


(22) 
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Figure  4.  Surface  before  (left)  and  after  (right)  low-pass  filter  at  the  moment  of  time 
t  ~  67To. 


where  D  =<  >  is  the  mean  square  amplitude. 

To  check  the  equation  (|22ll.  we  need  to  hnd  the  way  for  calculation  of  D{k).  If  the 
ensemble  is  stationary  in  time,  D{k)  could  be  found  for  any  given  k  by  the  time  averaging. 
In  our  case,  the  process  is  non-stationary,  and  we  have  a  problem  with  determination  of 
D{k). 

To  resolve  this  problem,  we  used  low-pass  hltering  instead  of  time  averaging.  The 
low-pass  hlter  was  taken  in  the  form 

/(^)  =  .  . 
A  =  0.25  ■  N^/2,  =  4096.  ^  ’ 

where  n  is  integer  vector  running  iF-space. 

This  choice  of  the  low-pass  hlter  preserves  the  values  of  total  energy,  wave  action  and  the 
total  momentum  within  three  percent  accuracy.  The  shape  of  low-pass  hltered  function 
is  presented  on  FigJU 

Now  it  is  possible  to  average  the  PDF  over  different  areas  in  /c-space.  The  results  for 
two  different  moments  of  time  t  ~  67.1To  and  t  ~  925To  are  presented  in  Figj5]and  Figj6l 

The  thin  line  gives  PDF  after  averaging  over  dissipation  region  harmonics,  while  bold 
line  presents  averaging  over  the  non-dissipative  area  \k\  <  k^  =  1024.  One  can  see  that 
statistics  in  the  last  case  is  quite  close  to  the  Gaussian,  while  in  the  dissipation  region 
it  deviates  from  Gaussian  distribution.  However,  deviation  from  the  Guassianity  in  the 
dissipation  region  does  not  create  any  problems,  since  the  “dissipative”  harmonics  do  not 
contain  any  essential  amount  of  the  total  energy,  wave  action  and  momentum. 

One  should  remember,  that  the  bold  lines  in  the  Figj5]  and  Figl6]  are  the  results  of 
averaging  over  a  million  of  harmonics.  Among  them  there  is  a  population  of  “selected 
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Figure  5.  Probability  distribution  function  for  relative  squared  amplitudes  |afcp/  < 
lofcp  >.  f  ~  67To.  Thin  and  bold  lines  -  averaging  over  dissipation  and  dissipation- 
free  regions  of  iC-space  correspondingly. 


few”,  or  “oligarchs”,  with  the  amplitudes  exceeding  the  average  value  by  the  factor  of 
more  than  ten  times.  The  oligarchs  exist  because  our  grid  is  still  not  hne  enough.  The 
contribution  of  such  oligarchs  in  our  case  to  the  total  wave  action  does  not  exceed  4%. 
Fifteen  leading  oligarchs  at  some  moment  of  time  are  presented  in  the  Appendix 
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Figure  6.  Probability  distribution  function  for  relative  squared  amplitudes  la^p/  < 
lofcp  >.  f  ~  925To.  Thin  and  bold  lines  -  averaging  over  dissipation  and  dissipation- 
free  regions  of  -ft'-space  correspondingly. 
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3.5.  Two-stage  evolution  of  the  swell 

Fig.  ITMTOl  demonstrate  time  evolution  of  main  characteristics  of  the  wave  field:  wave 
action,  energy,  characteristic  slope  and  mean  frequency. 


0 


Initial  Wave  Periods 

1000  2000  3000 


Solid  line  -  dynamical,  dashed  -  WAM3,  dotted  -  WAM4,  dash-dotted  -  artificial 


Figure  7.  Total  wave  action  as  a  function  of  time.  Solid  line  —  dynamical  equations, 
dashed-dotted  line  —  kinetic  equation  with  artificial  viscosity,  dashed  line  —  kinetic 
equation  with  WAM3  damping  term,  dotted  line  —  kinetic  equation  with  WAM4  damping 
term 


Figj9]  should  be  specially  commented.  Here  and  further  we  define  the  characteristic 
slope  as  defined  in  Eq.  HT]  According  to  this  definition  of  steepness  for  the  classical 
Pierson-Moscowitz  spectrum  /r  =  0.095.  Our  initial  steepness  /i  ~  0.15  exceeds  this  value 
essentially. 

Observed  evolution  of  the  spectrum  can  be  conventionally  separated  in  two  phases.  On 
the  first  stage  we  observe  fast  drop  of  wave  action,  slope  and  especially  energy.  Then 
the  drop  is  stabilized,  and  we  observe  slow  down-shift  of  mean  frequency  together  with 
angular  spreading.  Level  lines  of  smoothed  spectra  in  the  first  and  in  the  last  moments 
of  time  are  shown  in  FigJTTlIT^ 

Presence  of  two  stages  can  be  understood  through  the  study  of  the  Probability  Distri¬ 
bution  Functions  (PDFs)  of  the  elevations  of  the  surface.  In  all  figures  we  compare  the 
distribution  in  experimental  results  with  Gaussian  distribution  and  Tayfun  distribution  [ 
133].  The  later  case  is  just  a  first  correction  to  Gaussian  distribution  due  to  small  nonlin¬ 
earity.  We  used  explicit  form  of  Tayfun  distribution  following  [US]-  In  the  initial  moment 
of  time  the  PDF  is  Gaussian  (see  FigJT3|).  No  nonlinear  interaction  is  involved,  so  Tayfun 
distribution  does  not  fit  at  all.  However,  very  soon  intensive  super-Gaussian  tails  appear 
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Figure  8.  Total  wave  energy  as  a  function  of  time.  Solid  line  —  dynamical  equations, 
dashed-dotted  line  —  kinetic  equation  with  artificial  viscosity,  dashed  line  —  kinetic 
equation  with  WAM3  damping  term,  dotted  line  —  kinetic  equation  with  WAM4  damping 
term 


(see  FigJTTj).  They  are  well  described  by  Tayfun  distribution.  Then  tails  decrease  slowly 
(FigJT5|l.  and  in  the  last  moment  of  simulation,  when  characteristics  of  the  sea  are  close 
to  Peirson-Moscowitz,  the  statistics  is  close  to  Gaussian  again  (FigJT6|l.  Moderate  tails  do 
exist  and,  in  a  good  agreement  with  Tayfun  correction,  troughs  are  more  probable  than 
crests  which  in  turn  can  be  much  larger  in  absolute  value.  PDF  for  T]y  —  longitudinal 
gradients  in  the  first  moments  of  time  is  Gaussian  fFigJT7|).  Then  in  a  very  short  period 
of  time  strong  non-Gaussian  tails  appear  and  reach  their  maximum  at  f  ~  14To  (FigJTSl). 
Here  Tq  =  2ti/ —  period  of  initial  leading  wave.  Since  this  moment  the  non-Gaussian 
tails  decrease.  In  the  last  moment  of  simulation  they  are  essentially  reduced(FigJT^. 

Fast  growth  of  non-Gaussian  tails  can  be  explained  by  formation  of  coherent  harmonics. 
Indeed,  14To  ~  27r/(a;o/r)  is  a  characteristic  time  of  nonlinear  processes  due  to  quadratic 
nonlinearity.  Note  that  the  fourth  harmonic  in  our  system  is  strongly  decaying,  hence  we 
cannot  see  real  white  caps. 

Figures  [2(1112^  present  PDFs  for  gradients  in  the  orthogonal  direction. 
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Figure  9.  Average  wave  slope  as  a  function  of  time.  Solid  line  —  dynamical  equations, 
dashed-dotted  line  —  kinetic  equation  with  artihcial  viscosity,  dashed  line  —  kinetic 
eqnation  with  WAM3  damping  term,  dotted  line  —  kinetic  equation  with  WAM4  damping 
term 
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Solid  line  -  dynamical,  dashed  -  WAM3,  dotted  -  WAM4,  dash-dotted  -  artificial 


Figure  10.  Mean  wave  frequency  as  a  function  of  time.  Solid  line  —  dynamical  equations, 
dashed-dotted  line  —  kinetic  equation  with  artificial  viscosity,  dashed  line  —  kinetic 
equation  with  WAM3  damping  term,  dotted  line  —  kinetic  eqnation  with  WAM4  damping 
term 
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Figure  12.  Final  spectrum  |a^p.  t  =  3378To. 


18 


A.  O.  Korotkevich,  A.  Pushkarev,  D.  Resio  and  V.  E.  Zakharov 


Figure  13.  PDF  for  the  surface  elevation  rj  at  the  initial  moment  of  time,  t  =  0. 


Figure  14.  PDF  for  the  surface  elevation  rj  at  the  moment  of  maximum  surface  roughness. 
t  ~  14To. 


Numerical  Verification  of  the  Weak  Turbulent  Model  for  Swell  Evolution. 


19 


>* 


CC 

O 

D) 

O 

_l 


Figure  15.  PDF  for  the  surface  elevation  rj  at  some  middle  moment  of  time,  t  ~  67.1Tq. 
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Figure  16.  PDF  for  the  surface  elevation  r]  at  the  final  moment  of  time,  t  =  3378To. 
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Figure  21.  PDF  for  {'Vr])x  at  the  moment  of  maximum  surface  roughness,  t  ~  14To. 


Figure  22.  PDF  for  (Vri)^  at  the  hnal  moment  of  time,  t  =  3378To. 
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Figures  123112^  present  snapshots  of  the  surface  in  the  initial  and  hnal  moments  of 
simulation.  Figl25]  is  a  snapshot  of  the  surface  in  the  moment  of  maximal  roughness 
T  =  4.94  ~  14To. 
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Figure  23.  Surface  elevation  at  the  initial  moment  of  time,  t  =  0. 

4.  Statistical  numerical  experiment 

4.1.  Numerical  model  for  Hasselmann  Equation 

Numerical  integration  of  kinetic  equation  for  gravity  waves  on  deep  water  (Hasselmann 


equation)  was  the  subject  of  considerable  efforts  for  last  three  decades.  The  “ultimate 


goal”  of  the  effort  -  creation  of  the  operational  wave  model  for  wave  forecast  based 
on  direct  solution  of  the  Hasselmann  equation  -  happened  to  be  an  extremely  difficult 
computational  problem  due  to  mathematical  complexity  of  the  Sni  term,  which  requires 
calculation  of  a  three-dimensional  integral  at  every  moment  of  time. 

Historically,  numerical  methods  of  integration  of  kinetic  equation  for  gravity  waves  exist 
in  two  “flavors” . 

The  hrst  one  is  associated  with  works  [EU,  [135],  [EE],  [E7],  m  and  [  139] ,  and  is  based 
on  transformation  of  6-fold  into  3-fold  integrals  using  ^-functions.  Such  transformation 
leads  to  appearance  of  integrable  singularities,  which  creates  additional  difficulties  in 
calculations  of  the  Sni  term. 

The  second  type  of  models  has  been  developed  in  works  of  [  00]  and  [  01] ,  [  02]  and  is 
currently  known  as  Resio- Tracy  model.  It  uses  direct  calculation  of  resonant  quadruplet 
contribution  into  Sni  integral,  based  on  the  following  property:  given  two  hxed  vectors 
k,ki,  another  two  ^2,  ^3  are  uniquely  dehned  by  the  point  “moving”  along  the  resonant 
curve  -  locus. 
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Figure  24.  Surface  elevation  at  the  final  moment  of  time,  t  =  3378To. 

Numerical  simulation  in  the  current  work  was  performed  with  the  help  of  modihed 
version  of  the  second  type  algorithm.  Calculations  were  made  on  the  grid  71  x  36  points  in 


the  frequency-angle  domain  [cu,  6]  with  exponential  distribution  of  points  in  the  frequency 


domain  and  uniform  distribution  of  points  in  the  angle  direction. 

To  date,  Resio- Tracy  model  suffered  rigorous  testing  and  is  currently  used  with  high 
degree  of  trustworthiness:  it  was  tested  with  respect  to  motion  integrals  conservation  in 
the  “clean”  tests,  wave  action  conservation  in  wave  spectrum  down-shift,  realization  of 
self  -  similar  solution  in  “pure  swell”  and  “wind  forced”  regimes  (see  [SU,  [S3],  [S5]). 

Description  of  scaling  procedure  from  dynamical  equations  to  Hasselman  kinetic  equa¬ 
tion  variables  is  presented  in  Appendix  [Bl 

4.2.  Statistical  model  setup 

The  numerical  model  used  for  solution  of  the  Hasselmann  equation  has  been  supplied 
with  the  damping  term  in  three  different  forms: 

1.  Pseudo- viscous  high  frequency  damping  flTTll  used  in  dynamical  equations; 

2.  WAM3  viscous  term; 

3.  WAM4  viscous  term; 

Two  last  viscous  terms  are  the  “white-capping”  terms,  describing  energy  dissipation  by 
surface  waves  due  to  white-capping,  as  used  in  W AM  wave  forecasting  models  (in  SWAN 
model  the  used  term  has  different  but  close  tunable  parameters),  see  [146]: 


V 


(24) 


where  k  and  uj  are  wave  number  and  frequency,  tilde  denotes  mean  value;  Cds-,  <5  and  p 


are  tunable  coefficients;  S  =  ky/H  is  the  overall  steepness;  SpM  =  (3.02  x  10  is 
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Figure  25.  Surface  elevation  at  the  moment  of  maximum  surface  roughness,  t  ~  14To. 
Gradients  are  more  conspicuous. 


the  value  of  S  for  the  Pierson- Moscowitz  spectrum  (note  that  the  characteristic  steepness 
=  V2S).  It  is  worth  to  note  that  according  to  [|5T]  theoretical  value  of  steepness  for 
Pierson-Moscovitz  spectrum  is  SpM  —  (4.57  x  10“^)^'^.  It  gives  us  /r  ~  0.095. 

Values  of  tunable  coefficients  for  WAM3  case  are: 


Cds  =  2.36  X  10■^  5  =  0,  p  =  4 

(25) 

and  for  WAM4  case  are: 

Grf,  =  4.10  X  10-^  5  =  0.5,  p  =  4 

(26) 

In  all  three  cases  we  used  as  initial  data  smoothed  (hltered)  spectra  (see  FigJl])  obtained 
in  the  dynamical  run  at  the  time  =  3.65mm  =  24.3  ~  67.1To,  which  can  be  considered 
as  a  moment  of  the  end  of  the  hrst  “fast”  stage  of  spectral  evolution. 

The  natural  question  stemming  in  this  point,  is  why  calculation  of  the  dynamical  and 
Hasselmann  model  cannot  be  started  from  the  initial  conditions  flT6|)  simultaneously? 

There  are  good  reasons  for  that: 

As  it  was  mentioned  before,  the  time  evolution  of  the  initial  conditions  (1T6|1  in  presence 
of  the  damping  term  can  be  separated  in  two  stages:  relatively  fast  total  energy  drop 
in  the  beginning  of  the  evolution  and  succeeding  relatively  slow  total  energy  decrease  as 
a  function  of  time,  see  FigJHl  We  explain  this  phenomenon  by  existence  of  the  effective 
channel  of  the  energy  dissipation  due  to  strong  nonlinear  effects,  which  can  be  associated 
with  the  white-capping  as  we  mentioned  in  the  introduction. 

We  have  started  with  relatively  steep  waves  fi  ~  0.167.  As  we  see,  at  that  steepness 
white-capping  is  the  leading  effect.  This  fact  is  conhrmed  by  numerous  held  and  labora¬ 
tory  experiments.  From  the  mathematical  view-point,  the  white-capping  is  formation  of 
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coherent  structures  -  strongly  correlated  multiple  harmonics.  The  spectral  peak  is  posed 
in  our  experiments  initially  at  fc  ~  300,  while  the  edge  of  the  damping  area  kd  —  1024. 
Hence,  only  the  second  and  the  third  harmonic  can  be  developed,  while  higher  harmonics 
are  suppressed  by  strong  dissipation.  Anyway,  even  formation  of  the  second  and  the  third 
harmonic  is  enough  to  create  intensive  non-Gaussian  tail  of  the  PDF  for  longitudinal 
gradients.  This  process  is  very  fast.  In  the  moment  of  time  T  =  14To  we  see  fully  de¬ 
veloped  tails.  Relatively  sharp  gradients  mimic  formation  of  white  caps.  Certainly,  the 
“pure”  Hasselmann  equation  is  not  applicable  on  this  early  stage  of  spectral  evolution, 
when  energy  intensively  dissipates. 

As  steepness  decreases  and  spectral  maximum  of  the  swell  down-shifts,  the  efficiency  of 
such  mechanism  of  energy  absorption  becomes  less  important.  When  the  steepness  value 
drops  down  to  /i  ~  0.1,  at  approximately  T  ~  280To,  the  white-capping  is  negligibly 
small.  Therefore,  we  decided  to  start  comparison  between  deterministic  and  statistical 
modeling  in  some  intermediate  moment  of  time  T  ~  67.1To. 


5.  Comparison  of  deterministic  and  statistical  experiments. 

5.1.  Statistical  experiment  with  pseudo-viscous  damping  term. 

The  hrst  series  of  statistical  experiments  has  been  performed  with  pseudo- viscous  damp¬ 
ing  term  flTTl). 

FigJT]  -  [lO]  show  total  wave  action,  total  energy,  mean  wave  slope  and  mean  wave 
frequency  as  the  functions  of  time. 

FigJST]  shows  the  time  evolution  of  angle-averaged  wave  action  spectra  as  the  functions 
of  frequency  for  dynamical  and  Hasselmann  equations.  We  see  similar  down-shift  of  the 
spectral  maximum  both  in  dynamic  and  Hasselmann  equations.  The  correspondence  of 
the  spectral  maxima  amplitudes  is  not  good  at  all. 

It  is  quite  obvious  that  the  influence  of  the  artihcial  viscosity  is  not  strong  enough  to 
reach  the  correspondence  of  two  models. 

5.2.  Statistical  experiments  with  WAM3  damping  term 

The  second  series  of  statistical  experiments  has  been  done  for  the  choice  of  WAM3 
damping  term. 

The  temporal  behavior  of  total  wave  action,  energy  and  average  wave  slope  (see  FigJ3 
-  [TOl)  for  WAM3  damping  term  is  in  better  correspondence  with  dynamical  model,  than 
in  the  case  of  artihcial  viscosity  term.  For  initial  50  mm  duration  of  the  experiment  we 
observe  decent  correspondence  between  dynamical  and  Hasselmann  equations.  For  longer 
time  the  WAM3  model,  however,  deviates  from  the  dynamical  model  signihcantly. 

As  in  the  artihcial  viscosity  case,  the  angle-averaged  wave  action  spectra  as  the  func¬ 
tion  of  frequency  exhibit  distinct  down-shift  of  the  spectral  maxima  for  both  dynamical 
and  Hasselmann  equations  (see  FigJ32|).  Correspondence  between  time  evolution  of  the 
amplitudes  of  the  spectral  maxima  is  also  much  better  for  WAM3  choice  of  damping,  than 
for  artihcial  viscosity  case. 

Presumably,  WAM3  damping  term  underestimate  the  ehects  of  real  damping  at  the  very 
beginning  of  the  evolution  (when  the  ehects  of  white  capping  are  relatively  important), 
and  overestimates  them  on  later  stages  of  swell  evolution. 
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5.3.  Statistical  experiments  with  WAM4  damping  term 

The  final  third  series  of  experiments  have  been  done  for  the  choice  of  WAM4  damping 
term. 

FigJTHTO]  show  temporal  evolntion  of  the  total  wave  action,  total  energy,  mean  wave 
slope  and  mean  wave  frequency,  which  are  divergent  in  time  in  this  case. 

Figj33]  show  time  evolution  of  angle-averaged  wave  action  spectra  as  the  functions  of 
frequency  for  dynamical  and  Hasselmann  equations.  While  as  in  the  artificial  viscosity 
and  WAM3  cases  we  also  observe  distinct  down-shift  of  the  spectral  maxima,  the  corre¬ 
spondence  of  the  time  evolution  of  the  amplitudes  of  the  spectral  maxima  is  worse  than 
in  WAM3  case. 

This  observation  is  especially  surprising  in  view  of  the  fact  that  historically  WAM4 
damping  has  been  invented  as  an  improvement  to  WAM3  damping  term.  It  is  quite 
obvious  that  WAM4  damping  is  too  strong  for  description  of  the  reality  at  all  stages  of 
the  swell  evolution. 

6.  Down-shift  and  angnlar  spreading 

The  major  process  of  time-evolution  of  the  swell  is  frequency  down-shift.  During  T  = 
933To  the  mean  frequency  has  been  decreased  from  ujq  =  2  to  uji  =  0.6.  On  the  last  stage 
of  the  process,  the  mean  frequency  slowly  decays  as 

<  0^  >~  ~  t-1/15  (27) 

The  Hasselmann  equation  has  self-similar  solution,  describing  the  evolution  of  the  swell 
n{k,t)  ~  t^A^p  (^zTu)  (see  [03],  [05]).  For  this  solution 

<  0^  >~  t-A^^  (28) 

The  difference  between  fl^  and  (|28ll  can  be  explained  as  follows.  What  we  observed,  is 
not  a  self-similar  behavior.  Indeed,  a  self-similarity  presumes  that  the  angular  structure 
of  the  solution  is  constant  in  time.  Meanwhile,  we  observed  intensive  angular  spreading 
of  the  initially  narrow  in  angle,  almost  one- dimensional  wave  spectrum. 

Level  lines  of  the  low-pass  filtered  dynamic  and  kinetic  spectrum  at  the  beginning  of 
the  simulation  are  presented  on  FigJ26ll26l  There  are  ten  isolines  on  every  figure.  Values 
of  level  at  maximum  and  minimum  isolines  are  shown  on  every  picture.  Level  lines  of 
the  low-pass  filtered  dynamic  and  kinetic  spectrum  for  three  different  damping  terms  at 
the  end  of  the  simulation  are  presented  on  Figl?71l28l  We  observed  development  of 
bimodality  in  both  experiments.  This  is  in  accordance  with  field  observations  [|53],  [151]. 

One  can  see  good  correspondence  between  the  results  of  both  experiments.  Comparison 
of  time-evolution  of  the  mean  angular  spreading,  calculated  from  action  and  energy  spectra 
is  presented  on  Fig.  I29ll30l 
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Calculation  from  wave  action  spectrum 


Time(initial  wave  periods) 


Figure  29.  Average  angular  spreading  (^J  \d\n  {k)dk)  /  [j  n{k)dk)  as  a  function  of  time, 
calculated  from  wave  action.  Solid  line  -  WAM3,  dotted  line  -  WAM4,  dashed  line  - 
artificial  viscosity,  stars  -  dynamical  equations 


Calculation  from  wave  energy  spectrum 


Time(initial  wave  periods) 


Figure  30.  Average  angular  spreading  (^J  \6\un  (k)dk)  /  (/  (jn{k)d}Pj  as  a  function  of  time, 
calculated  from  wave  energy.  Solid  line  -  WAM3,  dotted  line  -  WAM4,  dashed  line  - 
artihcial  viscosity,  stars  -  dynamical  equations 
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We  see  growing  divergence  between  dynamic  and  kinetic  models.  But  using  WAM3  and 
WAM4  models  leads  to  worse  divergence.  This  is  an  additional  argument  against  these 
variants  of  white-cap  damping. 

One  can  expect  that  the  angular  spreading  will  be  arrested  at  later  times,  and  the 
spectra  will  take  universal  self-similar  shape. 


Figure  31.  Angle-averaged  spectrum  as  a  function  of  time  for  dynamical  and  Hasselmann 
equations  for  artihcial  viscosity  case. 


7.  Conclusion 

1.  Our  numerical  experiment  shows  that  the  Hasselmann  equation  without  any  ad¬ 
ditional  terms  is  applicable  for  description  of  gravity  waves  turbulence  in  the  in- 
hnite  basin  only  if  nonlinearity  is  low  enough.  The  average  wave  steepness  /i  = 
{E‘^ /N‘^)sqrt2E  should  be  signihcantly  less  a  certain  critical  level  /xq  =  0.095.  To 
describe  the  wave  turbulence  of  a  moderate  steepness  0.095  <  ja  <  0.15,  one  should 
add  to  the  Hasselmann  equation  some  additional  term  Sdiss  modeling  dissipation 
due  to  white-capping.  So  far  there  is  no  any  theory  making  possible  to  hnd  Sdiss 
analytically.  This  is  a  great  challenge  for  hydrodynamicists.  So  far  we  can  suggest 
that  Sdiss  depends  on  steepness  very  sharply.  Hence  the  guess  that  white-capping  is 
the  threshold  phenomenon  formulated  by  Banner,  Babanin,  and  Young  [[52]  looks 
very  plausible.  We  did  not  model  too  steep  waves,  but  we  can  conjecture  that  if 
/i  >  0.15,  white  capping  is  the  leading  nonlinear  process  and  the  weak-turbulent 
approach  is  not  applicable  at  ah 
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Figure  32.  Angle-averaged  spectrum  as  a  function  of  time  for  dynamical  and  Hasselmann 
equations  a  function  of  time  for  WAM3  case. 


2.  Our  results  are  not  a  surprise  for  designers  of  operational  models  for  wave  forecast¬ 
ing.  They  understood  necessity  to  introduce  into  the  models  extra  dissipative  terms 
many  years  ago.  Another  story  that  the  models  of  Sdiss  routinely  used  in  the  WAM 
and  WAVEWATCH  models  are  too  crude  and  do  not  grasp  the  most  important 
feature  of  the  white-capping  —  its  threshold  nature.  In  our  opinion  existing  models 
of  Sdiss  overestimate  dissipation  at  small  values  of  steepness  and  underestimate  it 
at  large  /i.  Moreover,  they  overestimate  dissipation  in  the  area  of  the  spectral  peak 
and  underestimate  it  in  the  spectral  region  of  high  wave  numbers.  We  are  planning 
to  offer  our  own  form  of  Sdiss  extracted  from  our  massive  numerical  simulation  of 
Euler  equation,  but  it  will  take  some  time  and  toil. 

3.  We  have  to  stress  that  we  solved  maximally  idealized  model.  We  did  not  take 
into  account  interaction  of  swell  with  atmosphere  (see  for  instance  [  |55],  [  |56]  ), 
interaction  with  ocean  currents,  deviation  of  the  wave  motion  from  potentiality  as 
well  as  influence  of  turbulence  in  the  surface-atmosphere  boundary  layer.  In  the 
future  we  plan  to  establish  a  contact  with  experimentalists  to  develop  more  realistic 
model  of  swell  propagation  in  the  real  ocean. 

4.  In  one  aspect  our  conclusions  are  very  resolute.  All  existing  experimental  wave  tanks 
cannot  be  used  for  modelling  of  the  waves  propagation  in  open  sea.  The  mesoscopic 
effects  in  wave  tanks  are  too  strong  for  reasonable  values  of  steepness.  This  pertains 
only  to  modeling  of  kinetics  of  energy  containing  region  in  the  vicinity  of  spectral 
peak.  The  rear  faces  of  spectral  distributions,  spectral  tails,  can  be  successfully 


34 


A.  O.  Korotkevich,  A.  Pushkarev,  D.  Resio  and  V.  E.  Zakharov 


Figure  33.  Angle-averaged  spectrum  as  a  function  of  time  for  dynamical  and  Hasselmann 
equations  a  function  of  time  for  WAM4  case. 


modelled.  This  was  demonstrated  by  Toba  in  his  classical  work  [|29].  Toba  observed 
the  Zakharov- Filonenko  spectrum  R  ~  in  the  wave  tank  of  moderate  size.  In 
our  opinion  this  conclusion  is  very  important.  Some  authors,  trying  to  hnd  the 
universal  law  of  fetch  dependence  of  mean  energy  and  mean  frequency  put  together 
data  collected  in  ocean  and  in  the  experimental  tanks.  This  mixed  data  hardly  can 
be  compared  with  any  self  sufficient  theory.  This  question  is  discussed  in  details  in 
the  article  \w\. 

Another  conclusion  is  more  pessimistic.  The  results  of  numerical  experiments  show  that  it 
is  very  difficult  to  reproduce  real  ocean  conditions  in  laboratory  wave  tank.  Even  the  size 
of  the  tanks  of  200  X  200  meters  is  not  large  enough  to  model  ocean  due  to  the  presence 
of  wave  numbers  grid  discreteness. 
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A.  “Forbes”  list  of  15  largest  harmonics. 

Here  one  can  find  15  largest  harmonics  at  the  end  of  calculations  in  the  framework  of 
dynamical  equations.  Average  square  of  amplitudes  in  dissipation-less  region  was  taken 
from  smoothed  spectrum,  while  all  these  harmonics  exceed  level  =  1.4  x  10“^^. 
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B.  Prom  Dynamical  Equations  to  Hasselmann  Equation. 


Standard  setup  for  numerical  simulation  of  the  dynamical  equations  ([1])  implies  27r  x  27r 
domain  in  real  space  and  gravity  acceleration  g  =  1.  Usage  of  the  domain  size  equal  27r 
is  convenient  because  in  this  case  wave  numbers  are  integers. 

In  the  contrary  to  dynamical  equations,  the  kinetic  equation  (IT^  is  formulated  in 
terms  of  real  physical  variables  and  it  is  necessary  to  describe  the  transformation  from 
the  “dynamical”  variables  into  to  the  “physical”  ones. 

EqJH  are  invariant  with  respect  to  “stretching”  transformation  from  “dynamical”  to 
“real”  variables: 


Vr 


t 


^  Ti  -> 

—k  ,  r  =  ar  , 

a 


9  =  vg', 


Ly  —  ah' 


(29) 

(30) 


where  prime  denotes  variables  corresponding  to  dynamical  equations. 

In  current  simulation  we  used  the  stretching  coefficient  a  =  800.00,  which  allows 
to  reformulate  the  statement  of  the  problem  in  terms  of  real  physics:  we  considered 
5026  m  X  5026  m  periodic  boundary  conditions  domain  of  statistically  uniform  ocean  with 
the  same  resolution  in  both  directions  and  characteristic  wave  length  of  the  initial  condi¬ 
tion  around  22  m.  In  oceanographic  terms,  this  statement  corresponds  to  the  “duration- 
limited  experiment” . 
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